rm(list = ls())
library(tidyr)
library(readr)
library(dplyr)
library(readxl)
library(tidyr)
library(arm)
library(lme4)
library(MASS)
library(texreg)
load("Data/paperdata.RData") #1993-2006

source("Code/coefplot.R")
summary(paperdata)

## DV: 
# ext_count
#ext_x_count
#ext_p
# ext_y
#ext_y_coun
#####
dv <- "ext_sup"
ctrs <- c("polity2_lag1","gdp_log_lag1", "military_log_lag1", "population_log_lag1",
          "bd_best_lag1", "deaths_civilians_lag1", "factor(gwno_loc)", "factor(year)")
## controls
ivs_1 <- c("troops_backward_un_lag1")
ivs_2 <- c("police_backward_un_lag1")
ivs_3 <- c("observers_backward_un_lag1")
ivs_4 <- c("totals_backward_un_lag1")


f_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
ivs <- c("troops_backward_un_lag1","police_backward_un_lag1","observers_backward_un_lag1","totals_backward_un_lag1")


## modeling
model_ext_sup_un <- list()
for (i in 1:length(ivs)){
  model_ext_sup_un[[i]] <- glm(eval(parse(text = paste0('f_', i))), 
                               data = paperdata, family = binomial(link = "logit"))
}

p2 <- marginaleffect_logit(ModelResults = model_ext_sup_un, n.sim = 1000, 
                           varname = c("troops_backward_un_lag1","police_backward_un_lag1",
                                       "observers_backward_un_lag1","totals_backward_un_lag1"),
                           data = paperdata,
                           val1 = 0, val2 =1, clusterid = "gwno_loc") 

p2 + scale_x_continuous(breaks = 1:4,
                        labels = parse(text =c( 
                          "维和部队",
                          "维和警察",
                          "维和观察员",
                          "所有类型")))+ 
  ggtitle("因变量：是否存在外部支持(联合国维和)")+
  labs(tag = "图3a ")+
  theme(plot.tag.position = "bottom")
# ggsave("Results/figures/fd_ext_sup_un.pdf", width = 6, height = 4)
# ggsave("Results/figures/fd_ext_sup_un.jpeg", dpi = 300, width = 8, height = 5, units = "cm")
ggsave("Replication/figures/figure3-a.bmp", width = 6, height = 4)
ggsave("Replication/figures/figure3-a.pdf", width = 6, height = 4)



### 支持政府或反政府


#####
dv <- "ext_sup_gov"
## controls
ivs_1 <- c("troops_backward_un_lag1")
ivs_2 <- c("police_backward_un_lag1")
ivs_3 <- c("observers_backward_un_lag1")
ivs_4 <- c("totals_backward_un_lag1")


f_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
ivs <- c("troops_backward_un_lag1","police_backward_un_lag1","observers_backward_un_lag1","totals_backward_un_lag1")


## modeling
model_ext_supgov_un <- list()
for (i in 1:length(ivs)){
  model_ext_supgov_un[[i]] <- glm(eval(parse(text = paste0('f_', i))), 
                                  data = paperdata, family = binomial(link = "logit"))
}

#####
dv <- "ext_sup_rebel"
## controls
ivs_1 <- c("troops_backward_un_lag1")
ivs_2 <- c("police_backward_un_lag1")
ivs_3 <- c("observers_backward_un_lag1")
ivs_4 <- c("totals_backward_un_lag1")


f_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
ivs <- c("troops_backward_un_lag1","police_backward_un_lag1","observers_backward_un_lag1","totals_backward_un_lag1")


## modeling
model_ext_suprebel_un <- list()
for (i in 1:length(ivs)){
  model_ext_suprebel_un[[i]] <- glm(eval(parse(text = paste0('f_', i))), 
                                    data = paperdata, family = binomial(link = "logit"))
}


p1_1 <- marginaleffect_logit2(ModelResults =  model_ext_supgov_un , n.sim = 1000, 
                              varname = c("troops_backward_un_lag1","police_backward_un_lag1",
                                          "observers_backward_un_lag1","totals_backward_un_lag1"),
                              data = paperdata,
                              dv = "ext_sup_gov",
                              val1 = 0, val2 =1, clusterid = "gwno_loc") 


p1_2 <- marginaleffect_logit2(ModelResults = model_ext_suprebel_un, n.sim = 1000, 
                              varname = c("troops_backward_un_lag1","police_backward_un_lag1",
                                          "observers_backward_un_lag1","totals_backward_un_lag1"),
                              data = paperdata,
                              dv = "ext_sup_rebel",
                              val1 = 0, val2 =1, clusterid = "gwno_loc") 




df_ext_sup_UN <- bind_rows(p1_1, p1_2)



p = ggplot(df_ext_sup_UN, aes(colour = dv)) + 
  geom_hline(yintercept = 0, lty = 2) +
  geom_linerange(aes(x = id, ymin = low,
                     ymax = high),
                 lwd = 1.5, position = position_dodge(width = 1/2)) +
  geom_pointrange(aes(x = id, y =median, ymin =low,
                      ymax =high, shape = dv),
                  fatten = 7,
                  lwd = 1.5, position = position_dodge(width = 1/2),
                  fill="white")+
  geom_text(aes(x = id +0.2 ,
                y = median,
                label = format(round(median, 2))),position = position_dodge(width = 1/2)) +
  theme_bw() + xlab("") + ylab("预测概率变化（边际效应）") + 
  theme(legend.position="bottom",
        legend.title=element_blank(),
        legend.text = element_text(colour="black", size=16),
        axis.text= element_text(size=16),
        text = element_text(size=16, family = 'KaiTi',face = "bold"),
        plot.title = element_text(hjust = .5, size = 16))+
  ggtitle("因变量：是否存在外部支持（联合国维和）")+
  labs(tag = "图3b ")+
  theme(plot.tag.position = "bottom")+
  scale_colour_manual(values = c("black", "black"), labels = c("支持政府","支持反叛组织"))+
  scale_shape_manual(values = c(15,23), labels = c("支持政府","支持反叛组织"))+
  guides(color = guide_legend(nrow = 1))+
  scale_x_continuous(breaks = 1:4,
                     labels = parse(text =c( 
                       "维和部队",
                       "维和警察",
                       "维和观察员",
                       "所有类型")))
#ggsave("Results/figures/fd_ext_sup_type_UN.pdf", width = 6, height = 4)
ggsave("Replication/figures/figure3-b.bmp", width = 6, height = 4)
ggsave("Replication/figures/figure3-b.pdf", width = 6, height = 4)


save(model_ext_sup_un,model_ext_supgov_un,model_ext_suprebel_un, file = "Results/Figure1nodels_UN.RData")



modSD = cluster_se(ModelResults = c(model_ext_sup_un,model_ext_supgov_un,model_ext_suprebel_un), data = paperdata, 
                   clusterid = "gwno_loc") 
modPvalue = cluster_pvalue(ModelResults = c(model_ext_sup_un,model_ext_supgov_un,model_ext_suprebel_un), data = paperdata, 
                           clusterid = "gwno_loc") 


screenreg(c(model_ext_sup_un,model_ext_supgov_un,model_ext_suprebel_un ))

htmlreg(c(model_ext_sup_un,model_ext_supgov_un,model_ext_suprebel_un ), file = "Replication/tables/ext_sup_UN.doc",
        stars = c(0.01, 0.05, 0.1),
        override.se =  modSD,
        override.pvalues = modPvalue,
        caption = "Logit回归结果：联合国维和行动对外部支持的影响",
        caption.above = TRUE,
        label = "tab:tab5",
        scalebox = 0.7,
        custom.header = list("外部支持" = 1:4,
                             "支持政府" = 5:8,
                             "支持反叛组织" = 9:12),
        # custom.model.names = c("M1:troop", "M2:police", "M3:observers", "M4:any",
        #                        "M5:troop", "M6:police", "M7:observers", "M8:any",
        #                        "M9:troop", "M10:police", "M11:observers", "M12:any"),
        custom.coef.map = list("troops_backward_un_lag1" = "维和部队",
                               "police_backward_un_lag1" = "维和警察",
                               "observers_backward_un_lag1" ="维和观察员",
                               "totals_backward_un_lag1" = "所有类型",
                               "polity2_lag1" = "政体类型",
                               "gdp_log_lag1" = "国民生产总值（对数）",
                               "military_log_lag1" = "军队人数（对数）",
                               "population_log_lag1" = "人口总数(对数)",
                               "bd_best_lag1"= "战斗相关死亡人数",
                               "deaths_civilians_lag1" = "平民死亡人数"),
        custom.gof.rows = list("国家固定效应" = c("是", "是", "是", "是",
                                            "是", "是", "是", "是",
                                            "是", "是", "是", "是"),
                               "年份固定效应" = c("是", "是", "是", "是",
                                            "是", "是", "是", "是",
                                            "是", "是", "是", "是")),
        digits = 3)

